Generating final model predictions after relativizing
Alice Stears
2025-11-12
##User defined parameters and functions
print(params)## $readParams
## [1] TRUE# set to true if want to run for a limited number of rows (i.e. for code testing)
readParams <- params$readParams
makePredictions <- function(predictionDF, modelObject) {
##
# predictionDF <- climDatPred
# modelObject <- bestLambdaMod_grassShrub_totalHerb
# ##
# pretend to scale the input variables so the model object can predict accurately
predictionDF <- predictionDF %>%
mutate(across(all_of(prednames), base::scale,scale = FALSE, center = FALSE))
# modelPredictions
modelPreds <- predict(object = modelObject, newdata = predictionDF, type = "response")
# add predictions back into the input data.frame
predictionDF <- predictionDF %>%
cbind(modelPreds)
# truncate all predictions to max out at 100
#predictionDF[predictionDF$modelPreds>100 & !is.na(predictionDF$modelPreds),"modelPreds"] <- 100
predictionDF[predictionDF$modelPreds < 0 & !is.na(predictionDF$modelPreds),"modelPreds"] <- 0
# print predicted data
return(predictionDF)
}Dependencies
library(tidyverse)
library(sf)
library(terra)
library(kableExtra)
library(knitr)
library(USA.state.boundaries)
library(tidyterra)
library(ggpubr)
library(USA.state.boundaries)Read in necessary data
# Load Data ---------------------------------------------------------------
# data ready for modeling (that has been scaled)
modDat_1_s <- readRDS("./models/scaledModelInputData.rds")
# get the soil raster, which we'll use for rasterizing the imagery
soilRastTemp <- readRDS("../../../Data_processed/SoilsRaster.rds") %>%
terra::unwrap()
# make a map of the predictions
test_rast <- rast("../../../Data_raw/dayMet/rawMonthlyData/orders/70e0da02b9d2d6e8faa8c97d211f3546/Daymet_Monthly_V4R1/data/daymet_v4_prcp_monttl_na_1980.tif") %>%
terra::aggregate(fact = 12, fun = "mean")
# download map info for visualization
data(state_boundaries_wgs84)
cropped_states <- suppressMessages(state_boundaries_wgs84 %>%
dplyr::filter(NAME!="Hawaii") %>%
dplyr::filter(NAME!="Alaska") %>%
dplyr::filter(NAME!="Puerto Rico") %>%
dplyr::filter(NAME!="American Samoa") %>%
dplyr::filter(NAME!="Guam") %>%
dplyr::filter(NAME!="Commonwealth of the Northern Mariana Islands") %>%
dplyr::filter(NAME!="United States Virgin Islands") %>%
dplyr::filter(NAME!= "U.S. Virgin Islands") %>%
sf::st_sf() %>%
sf::st_transform(sf::st_crs(test_rast)))
#sf::st_crop(sf::st_bbox(modDat_1_sf)+c(-1,-1,1,1))
## add ecoregion boundaries (for our ecoregion level model)
regions <- sf::st_read(dsn = "../../../Data_raw/Level2Ecoregions/", layer = "NA_CEC_Eco_Level2") ## Reading layer `NA_CEC_Eco_Level2' from data source
## `/Users/astears/Documents/Dropbox_static/Work/NAU_USGS_postdoc/cleanPED/PED_vegClimModels/Data_raw/Level2Ecoregions'
## using driver `ESRI Shapefile'
## Simple feature collection with 2261 features and 8 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -4334052 ymin: -3313739 xmax: 3324076 ymax: 4267265
## Projected CRS: Sphere_ARC_INFO_Lambert_Azimuthal_Equal_Arearegions <- regions %>%
st_transform(crs = st_crs(test_rast)) %>%
st_make_valid()
ecoregionLU <- data.frame("NA_L1NAME" = sort(unique(regions$NA_L1NAME)),
"newRegion" = c(NA, "Forest", "dryShrubGrass",
"dryShrubGrass", "Forest", "dryShrubGrass",
"dryShrubGrass", "Forest", "Forest",
"dryShrubGrass", "Forest", "Forest",
"Forest", "Forest", "dryShrubGrass",
NA
))
goodRegions <- regions %>%
left_join(ecoregionLU)
mapRegions <- goodRegions %>%
filter(!is.na(newRegion)) %>%
group_by(newRegion) %>%
summarise(geometry = sf::st_union(geometry)) %>%
ungroup() %>%
st_simplify(dTolerance = 1000)
## contemporary climate data
climDat_temp <- readRDS("/Users/astears/Documents/Dropbox_static/Work/NAU_USGS_postdoc/PED_vegClimModels/Data_processed/EcoRegion_climSoilData.rds")
climDat_timeTemp <- climDat_temp %>%
mutate(spatialID = paste0(round(x,4), "_", round(y,4))) %>%
drop_na()
spatID_lu <- data.frame("spatialID" = unique(climDat_timeTemp$spatialID),
"uniqueID" = 1:length(unique(climDat_timeTemp$spatialID)))
climDat_timeTemp <- climDat_timeTemp %>%
left_join(spatID_lu)
## because of the spatial joining we did, there are some "locations" that have data for less than the true time span we want (2010-2020), or multiple observations per year... remove those so we have the even representation over the time span we're looking for
set.seed("12011993")
goodTable <- data.frame(table(climDat_timeTemp$uniqueID))
goodIDs_temp <- as.numeric(goodTable[goodTable$Freq == 11,"Var1"])
#goodIDs <- unique(names(table(climDat_timeTemp$uniqueID))[(table(climDat_timeTemp$uniqueID)==11)])
# select 1000 random IDs
goodIDs <- goodIDs_temp[sample(x = 1:length(goodIDs_temp), size = 1000, replace = FALSE)]
## get climate data for 500 different points randomly chosen across CONUS for all years in the contemporary climate/weather/soils dataset
climDat_time <- climDat_timeTemp %>%
filter(uniqueID %in% goodIDs) %>%
dplyr::select(tmin_meanAnnAvg_CLIM:durationFrostFreeDays_meanAnnAvg_3yrAnom, NA_L1CODE,
NA_L1NAME, NA_L1KEY, year, newRegion, x, y, soilDepth:totalAvailableWaterHoldingCapacity, uniqueID, spatialID) %>%
rename("tmin" = tmin_meanAnnAvg_CLIM,
"tmax" = tmax_meanAnnAvg_CLIM, #1
"tmean" = tmean_meanAnnAvg_CLIM,
"prcp" = prcp_meanAnnTotal_CLIM,
"t_warm" = T_warmestMonth_meanAnnAvg_CLIM,
"t_cold" = T_coldestMonth_meanAnnAvg_CLIM,
"prcp_wet" = precip_wettestMonth_meanAnnAvg_CLIM,
"prcp_dry" = precip_driestMonth_meanAnnAvg_CLIM,
"prcp_seasonality" = precip_Seasonality_meanAnnAvg_CLIM, #2
"prcpTempCorr" = PrecipTempCorr_meanAnnAvg_CLIM, #3
"abvFreezingMonth" = aboveFreezing_month_meanAnnAvg_CLIM,
"isothermality" = isothermality_meanAnnAvg_CLIM, #4
"annWatDef" = annWaterDeficit_meanAnnAvg_CLIM,
"annWetDegDays" = annWetDegDays_meanAnnAvg_CLIM,
"VPD_mean" = annVPD_mean_meanAnnAvg_CLIM,
"VPD_max" = annVPD_max_meanAnnAvg_CLIM, #5
"VPD_min" = annVPD_min_meanAnnAvg_CLIM, #6
"VPD_max_95" = annVPD_max_95percentile_CLIM,
"annWatDef_95" = annWaterDeficit_95percentile_CLIM,
"annWetDegDays_5" = annWetDegDays_5percentile_CLIM,
"frostFreeDays_5" = durationFrostFreeDays_5percentile_CLIM,
"frostFreeDays" = durationFrostFreeDays_meanAnnAvg_CLIM,
"soilDepth" = soilDepth, #7
"clay" = surfaceClay_perc,
"sand" = avgSandPerc_acrossDepth, #8
"coarse" = avgCoarsePerc_acrossDepth, #9
"carbon" = avgOrganicCarbonPerc_0_3cm, #10
"AWHC" = totalAvailableWaterHoldingCapacity,
## anomaly variables
tmean_anom = tmean_meanAnnAvg_3yrAnom, #15
tmin_anom = tmin_meanAnnAvg_3yrAnom, #16
tmax_anom = tmax_meanAnnAvg_3yrAnom, #17
prcp_anom = prcp_meanAnnTotal_3yrAnom, #18
t_warm_anom = T_warmestMonth_meanAnnAvg_3yrAnom, #19
t_cold_anom = T_coldestMonth_meanAnnAvg_3yrAnom, #20
prcp_wet_anom = precip_wettestMonth_meanAnnAvg_3yrAnom, #21
precp_dry_anom = precip_driestMonth_meanAnnAvg_3yrAnom, #22
prcp_seasonality_anom = precip_Seasonality_meanAnnAvg_3yrAnom, #23
prcpTempCorr_anom = PrecipTempCorr_meanAnnAvg_3yrAnom, #24
aboveFreezingMonth_anom = aboveFreezing_month_meanAnnAvg_3yrAnom, #25
isothermality_anom = isothermality_meanAnnAvg_3yrAnom, #26
annWatDef_anom = annWaterDeficit_meanAnnAvg_3yrAnom, #27
annWetDegDays_anom = annWetDegDays_meanAnnAvg_3yrAnom, #28
VPD_mean_anom = annVPD_mean_meanAnnAvg_3yrAnom, #29
VPD_min_anom = annVPD_min_meanAnnAvg_3yrAnom, #30
VPD_max_anom = annVPD_max_meanAnnAvg_3yrAnom, #31
VPD_max_95_anom = annVPD_max_95percentile_3yrAnom, #32
annWatDef_95_anom = annWaterDeficit_95percentile_3yrAnom, #33
annWetDegDays_5_anom = annWetDegDays_5percentile_3yrAnom , #34
frostFreeDays_5_anom = durationFrostFreeDays_5percentile_3yrAnom, #35
frostFreeDays_anom = durationFrostFreeDays_meanAnnAvg_3yrAnom #36
) %>%
dplyr::select(-c(tmin_meanAnnAvg_3yr:Start_3yr))
# rename
climDat <- climDat_temp %>%
filter(year == 2016) %>%
dplyr::select(tmin_meanAnnAvg_CLIM:durationFrostFreeDays_meanAnnAvg_3yrAnom, NA_L1CODE,
NA_L1NAME, NA_L1KEY, newRegion, x, y, soilDepth:totalAvailableWaterHoldingCapacity) %>%
rename("tmin" = tmin_meanAnnAvg_CLIM,
"tmax" = tmax_meanAnnAvg_CLIM, #1
"tmean" = tmean_meanAnnAvg_CLIM,
"prcp" = prcp_meanAnnTotal_CLIM,
"t_warm" = T_warmestMonth_meanAnnAvg_CLIM,
"t_cold" = T_coldestMonth_meanAnnAvg_CLIM,
"prcp_wet" = precip_wettestMonth_meanAnnAvg_CLIM,
"prcp_dry" = precip_driestMonth_meanAnnAvg_CLIM,
"prcp_seasonality" = precip_Seasonality_meanAnnAvg_CLIM, #2
"prcpTempCorr" = PrecipTempCorr_meanAnnAvg_CLIM, #3
"abvFreezingMonth" = aboveFreezing_month_meanAnnAvg_CLIM,
"isothermality" = isothermality_meanAnnAvg_CLIM, #4
"annWatDef" = annWaterDeficit_meanAnnAvg_CLIM,
"annWetDegDays" = annWetDegDays_meanAnnAvg_CLIM,
"VPD_mean" = annVPD_mean_meanAnnAvg_CLIM,
"VPD_max" = annVPD_max_meanAnnAvg_CLIM, #5
"VPD_min" = annVPD_min_meanAnnAvg_CLIM, #6
"VPD_max_95" = annVPD_max_95percentile_CLIM,
"annWatDef_95" = annWaterDeficit_95percentile_CLIM,
"annWetDegDays_5" = annWetDegDays_5percentile_CLIM,
"frostFreeDays_5" = durationFrostFreeDays_5percentile_CLIM,
"frostFreeDays" = durationFrostFreeDays_meanAnnAvg_CLIM,
"soilDepth" = soilDepth, #7
"clay" = surfaceClay_perc,
"sand" = avgSandPerc_acrossDepth, #8
"coarse" = avgCoarsePerc_acrossDepth, #9
"carbon" = avgOrganicCarbonPerc_0_3cm, #10
"AWHC" = totalAvailableWaterHoldingCapacity,
## anomaly variables
tmean_anom = tmean_meanAnnAvg_3yrAnom, #15
tmin_anom = tmin_meanAnnAvg_3yrAnom, #16
tmax_anom = tmax_meanAnnAvg_3yrAnom, #17
prcp_anom = prcp_meanAnnTotal_3yrAnom, #18
t_warm_anom = T_warmestMonth_meanAnnAvg_3yrAnom, #19
t_cold_anom = T_coldestMonth_meanAnnAvg_3yrAnom, #20
prcp_wet_anom = precip_wettestMonth_meanAnnAvg_3yrAnom, #21
precp_dry_anom = precip_driestMonth_meanAnnAvg_3yrAnom, #22
prcp_seasonality_anom = precip_Seasonality_meanAnnAvg_3yrAnom, #23
prcpTempCorr_anom = PrecipTempCorr_meanAnnAvg_3yrAnom, #24
aboveFreezingMonth_anom = aboveFreezing_month_meanAnnAvg_3yrAnom, #25
isothermality_anom = isothermality_meanAnnAvg_3yrAnom, #26
annWatDef_anom = annWaterDeficit_meanAnnAvg_3yrAnom, #27
annWetDegDays_anom = annWetDegDays_meanAnnAvg_3yrAnom, #28
VPD_mean_anom = annVPD_mean_meanAnnAvg_3yrAnom, #29
VPD_min_anom = annVPD_min_meanAnnAvg_3yrAnom, #30
VPD_max_anom = annVPD_max_meanAnnAvg_3yrAnom, #31
VPD_max_95_anom = annVPD_max_95percentile_3yrAnom, #32
annWatDef_95_anom = annWaterDeficit_95percentile_3yrAnom, #33
annWetDegDays_5_anom = annWetDegDays_5percentile_3yrAnom , #34
frostFreeDays_5_anom = durationFrostFreeDays_5percentile_3yrAnom, #35
frostFreeDays_anom = durationFrostFreeDays_meanAnnAvg_3yrAnom #36
) %>%
dplyr::select(-c(tmin_meanAnnAvg_3yr:Start_3yr))
rm(climDat_temp)
gc()## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 3027807 161.8 4932448 263.5 NA 3877882 207.2
## Vcells 950769060 7253.8 2664304804 20327.1 65536 2664034971 20325.0## now, scale the contemporary climate data for use in models
# get the scaling factors
scaleParams <- modDat_1_s %>%
#filter(Year == 2016) %>%
dplyr::select(tmin_s:AWHC_s) %>%
reframe(across(all_of(names(.)), attributes))
# apply the scaling factors to the contemporary climate data
namesToScale <- climDat %>%
dplyr::select(tmin:frostFreeDays, tmean_anom:frostFreeDays_anom, soilDepth:AWHC) %>%
names()
climDat_scaled <- map(namesToScale, .f = function(x) {
x_new <- (climDat[,x] - scaleParams[,paste0(x, "_s")]$`scaled:center`)/scaleParams[,paste0(x, "_s")]$`scaled:scale`
return(data.frame(x_new))
}) %>%
purrr::list_cbind()
names(climDat_scaled) <- paste0(namesToScale, "_s")
climDatPred <- climDat %>%
dplyr::select(NA_L1CODE:y) %>%
cbind(climDat_scaled)
names(climDatPred)[7:56] <- str_remove(names(climDatPred)[7:56], pattern = "_s$")
# apply the scaling factors to the contemporary climate data used for predictions over time
climDat_scaled <- map(namesToScale, .f = function(x) {
x_new <- (climDat_time[,x] - scaleParams[,paste0(x, "_s")]$`scaled:center`)/scaleParams[,paste0(x, "_s")]$`scaled:scale`
return(data.frame(x_new))
}) %>%
purrr::list_cbind()
names(climDat_scaled) <- paste0(namesToScale, "_s")
climDat_time_Pred <- climDat_time %>%
dplyr::select(NA_L1CODE:y,spatialID) %>%
cbind(climDat_scaled)
names(climDat_time_Pred)[9:58] <- str_remove(names(climDat_time_Pred)[9:58], pattern = "_s$")
rm(climDat_scaled)
gc()## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 3035195 162.1 4932448 263.5 NA 4932448 263.5
## Vcells 976334070 7448.9 2664304804 20327.1 65536 2664034971 20325.0prednames_s <- modDat_1_s %>%
dplyr::select(tmin_s:AWHC_s) %>%
names()
prednames <- str_replace(prednames_s, pattern = "_s$", replacement = "")
## read in scaled data for future climate
forecastClimSoilsDatPred_1 <- readRDS(file = "../../../Data_processed/CoverData/IntermediateAnalysisFiles/Final_ForecastedClimateDataAndSoilsDataForPredictions_BNU-ESM_rcp8_5.rds")
forecastClimSoilsDatPred_2 <- readRDS(file = "../../../Data_processed/CoverData/IntermediateAnalysisFiles/Final_ForecastedClimateDataAndSoilsDataForPredictions_IPSL-CM5A-MR_rcp8_5.rds")
# read in unscaled data for future climate
forecastClimSoilsDat_1 <- readRDS("../../../Data_processed/CoverData/IntermediateAnalysisFiles/Final_ForecastedClimateDataAndSoilsDataForPredictions_BNU-ESM_rcp8_5_UNSCALED.rds")
forecastClimSoilsDat_2 <- readRDS("../../../Data_processed/CoverData/IntermediateAnalysisFiles/Final_ForecastedClimateDataAndSoilsDataForPredictions_IPSL-CM5A-MR_rcp8_5_UNSCALED.rds")Now, predict using the entire hierarchy of models, starting with the ecoregion classification model
These are predictions of ecoregion classification using contemporary and modeled future climate data
ecoMod <- readRDS("../../ecoregionClassification/ModelFitting/ModelObjects/EcoregionClassificationModel.rds")
## print out model statement
coefficientTable <- data.frame(coefficients(ecoMod))
responseVar <- "P(forest)"
coefficientTable$coefficientName <- rownames(coefficientTable)
coefficientTable$coefficients.ecoMod. <- round(coefficientTable$coefficients.ecoMod., 4)
# print out coefficients w/ coefficient names
tempNames <- paste0(
apply(coefficientTable, MARGIN = 1, FUN = function(x) {
if (x["coefficientName"] == "(Intercept)") {
paste0(x["coefficients.ecoMod."])
} else {
paste0(x["coefficients.ecoMod."], "*", x["coefficientName"])
}
}
),
collapse = " + "
)
# print the unscaled model statement
unscaledModelName <- paste0(responseVar, "~ 1/ (1 + exp(-(", tempNames,")))")
print(unscaledModelName)## [1] "P(forest)~ 1/ (1 + exp(-(10.1438 + -0.3121*T_warmestMonth_meanAnnAvg_CLIM + 0.2527*T_coldestMonth_meanAnnAvg_CLIM + 0.0106*precip_wettestMonth_meanAnnAvg_CLIM + -0.0585*annWaterDeficit_meanAnnAvg_CLIM + -2.5967*PrecipTempCorr_meanAnnAvg_CLIM + 0.0582*isothermality_meanAnnAvg_CLIM + -0.0081*soilDepth + 0.0328*avgSandPerc_acrossDepth + 0.0304*avgCoarsePerc_acrossDepth + 0.2734*avgOrganicCarbonPerc_0_3cm)))"Here are maps of the ecoregion model predictions with both contemporary and modeled future climate data
## get model data used for ecoregion fitting
modDat_ecoregionFit <- readRDS("../../../Data_processed/EcoRegion_climSoilData.rds")
modDat_ecoregionFit <- modDat_ecoregionFit %>%
rename("Long" = x, "Lat" = y) %>%
dplyr::select(c(newRegion, #swe_meanAnnAvg_CLIM:
tmin_meanAnnAvg_CLIM:durationFrostFreeDays_meanAnnAvg_CLIM,
soilDepth , surfaceClay_perc ,
avgSandPerc_acrossDepth , avgCoarsePerc_acrossDepth,
avgOrganicCarbonPerc_0_3cm , totalAvailableWaterHoldingCapacity,
Long, Lat)) %>%
mutate(newRegion = as.factor(newRegion)) %>%
drop_na()
# get climate data from dayMet (d.f. is "climDatPred_unscaled")
climDatPred_unscaled <- climDat
names(climDatPred_unscaled)[c(5, 6, 7, 13, 10, 12, 98, 100, 101, 102)] <- c("T_warmestMonth_meanAnnAvg_CLIM", "T_coldestMonth_meanAnnAvg_CLIM", "precip_wettestMonth_meanAnnAvg_CLIM", "annWaterDeficit_meanAnnAvg_CLIM", "PrecipTempCorr_meanAnnAvg_CLIM", "isothermality_meanAnnAvg_CLIM", "soilDepth", "avgSandPerc_acrossDepth", "avgCoarsePerc_acrossDepth", "avgOrganicCarbonPerc_0_3cm")
#rm(climDat)
#gc()
# predict with contemporary climate data
preds_byHand <- climDatPred_unscaled %>%
mutate(pred = 1/(1 + exp(-( 9.8726 + -0.2999*T_warmestMonth_meanAnnAvg_CLIM + 0.2456*T_coldestMonth_meanAnnAvg_CLIM + 0.0106*precip_wettestMonth_meanAnnAvg_CLIM + -0.0621*annWaterDeficit_meanAnnAvg_CLIM + -2.7863*PrecipTempCorr_meanAnnAvg_CLIM + 0.0540*isothermality_meanAnnAvg_CLIM + -0.0076*soilDepth + 0.0335*avgSandPerc_acrossDepth + 0.0310*avgCoarsePerc_acrossDepth + 0.2726*avgOrganicCarbonPerc_0_3cm))))
#predictionsModel <- predict(object = ecoMod, type = "response", newdata = climDatPred_unscaled)
# predict with modeled future climate data v1
names(forecastClimSoilsDat_1)[c(8, 9, 10, 16, 13, 15, 48, 50, 51, 52)] <- c("T_warmestMonth_meanAnnAvg_CLIM", "T_coldestMonth_meanAnnAvg_CLIM", "precip_wettestMonth_meanAnnAvg_CLIM", "annWaterDeficit_meanAnnAvg_CLIM", "PrecipTempCorr_meanAnnAvg_CLIM", "isothermality_meanAnnAvg_CLIM", "soilDepth", "avgSandPerc_acrossDepth", "avgCoarsePerc_acrossDepth", "avgOrganicCarbonPerc_0_3cm")
preds_future1_byHand <- forecastClimSoilsDat_1 %>%
mutate(pred = 1/(1 + exp(-( 9.8726 + -0.2999*T_warmestMonth_meanAnnAvg_CLIM + 0.2456*T_coldestMonth_meanAnnAvg_CLIM + 0.0106*precip_wettestMonth_meanAnnAvg_CLIM + -0.0621*annWaterDeficit_meanAnnAvg_CLIM + -2.7863*PrecipTempCorr_meanAnnAvg_CLIM + 0.0540*isothermality_meanAnnAvg_CLIM + -0.0076*soilDepth + 0.0335*avgSandPerc_acrossDepth + 0.0310*avgCoarsePerc_acrossDepth + 0.2726*avgOrganicCarbonPerc_0_3cm))))
# predict with modeled future climate data v2
names(forecastClimSoilsDat_2)[c(8, 9, 10, 16, 13, 15, 48, 50, 51, 52)] <- c("T_warmestMonth_meanAnnAvg_CLIM", "T_coldestMonth_meanAnnAvg_CLIM", "precip_wettestMonth_meanAnnAvg_CLIM", "annWaterDeficit_meanAnnAvg_CLIM", "PrecipTempCorr_meanAnnAvg_CLIM", "isothermality_meanAnnAvg_CLIM", "soilDepth", "avgSandPerc_acrossDepth", "avgCoarsePerc_acrossDepth", "avgOrganicCarbonPerc_0_3cm")
preds_future2_byHand <- forecastClimSoilsDat_2 %>%
mutate(pred = 1/(1 + exp(-( 9.8726 + -0.2999*T_warmestMonth_meanAnnAvg_CLIM + 0.2456*T_coldestMonth_meanAnnAvg_CLIM + 0.0106*precip_wettestMonth_meanAnnAvg_CLIM + -0.0621*annWaterDeficit_meanAnnAvg_CLIM + -2.7863*PrecipTempCorr_meanAnnAvg_CLIM + 0.0540*isothermality_meanAnnAvg_CLIM + -0.0076*soilDepth + 0.0335*avgSandPerc_acrossDepth + 0.0310*avgCoarsePerc_acrossDepth + 0.2726*avgOrganicCarbonPerc_0_3cm))))
## Plot predictions with contemporary data
# make into a raster
plotObs <- preds_byHand %>%
#drop_na(paste(response)) %>%
#slice_sample(n = 5e4) %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = "pred",
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
# get the extent of this particular raster, and crop it accordingly
tempExt <- crds(plotObs, na.rm = TRUE)
plotObs_2 <- plotObs %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
map_ecoRegion_current <- ggplot() +
geom_spatraster(data = plotObs_2) +
geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Ecoregion classification predictions using
contemporary dayMet data")) +
scale_fill_gradient2(low = "brown",
mid = "wheat" ,
high = "darkgreen" ,
midpoint = 0, limits = c(0,1), na.value = "lightgrey",
labels = c("grassShrub", "0.25", "0.50", "0.75", "forest")) +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)])
# for future v1
plotObs_future1 <- preds_future1_byHand %>%
#drop_na(paste(response)) %>%
#slice_sample(n = 5e4) %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = "pred",
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
plotObs_future1_2 <- plotObs_future1 %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
map_ecoRegion_future1 <- ggplot() +
geom_spatraster(data = plotObs_future1_2) +
geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_future1_2)),fill=NA ) +
labs(title = paste0("Ecoregion classification predictions using future
climate data from the BNU-ESM model")) +
scale_fill_gradient2(low = "brown",
mid = "wheat" ,
high = "darkgreen" ,
midpoint = 0, limits = c(0,1), na.value = "lightgrey",
labels = c("grassShrub", "0.25", "0.50", "0.75", "forest")) +
xlim(st_bbox(plotObs_future1_2)[c(1,3)]) +
ylim(st_bbox(plotObs_future1_2)[c(2,4)])
# for future v2
plotObs_future2 <- preds_future2_byHand %>%
#drop_na(paste(response)) %>%
#slice_sample(n = 5e4) %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = "pred",
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
plotObs_future2_2 <- plotObs_future2 %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
map_ecoRegion_future2 <- ggplot() +
geom_spatraster(data = plotObs_future2_2) +
geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_future2_2)),fill=NA ) +
labs(title = paste0("Ecoregion classification predictions using future
climate data from the IPSL-CM5A-MR model")) +
scale_fill_gradient2(low = "brown",
mid = "wheat" ,
high = "darkgreen" ,
midpoint = 0, limits = c(0,1), na.value = "lightgrey",
labels = c("grassShrub", "0.25", "0.50", "0.75", "forest")) +
xlim(st_bbox(plotObs_future2_2)[c(1,3)]) +
ylim(st_bbox(plotObs_future2_2)[c(2,4)])
## Plot the maps together
ggarrange(map_ecoRegion_current, map_ecoRegion_future1, map_ecoRegion_future2, nrow = 1) %>%
annotate_figure(fig.lab = "Model Predictions of Ecoregion Classification with Contemporary and Forecasted Climate Data", fig.lab.size = 20)
get model objects for cover models
# Total herbaceous cover:
# Grass/shrub: 1SE lambda model
totalHerb_GS <- readRDS("./models/betaLASSO/TotalHerbaceousCover_shrubGrass_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_oneSELambdaGLM.rds")
# Forest: best lambda model
totalHerb_F <- readRDS("./models/betaLASSO/TotalHerbaceousCover_forest_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
# Total tree cover:
# Grass/shrub: best lambda model
totalTree_GS <- readRDS("./models/betaLASSO/TotalTreeCover_shrubGrass_noTLP_FALSE_removeAnomaliesTRUE_bestLambdaGLM.rds")
# Forest: best lambda model
totalTree_F <- readRDS("./models/betaLASSO/TotalTreeCover_forest_noTLP_FALSE_removeAnomaliesTRUE_halfSELambdaGLM.rds")
# Bare ground:
# CONUS - 1SE lambda model
bareGround_CONUS <- readRDS("./models/betaLASSO/BareGroundCover_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_oneSELambdaGLM.rds")
# Shrub:
# CONUS - best lambda model
shrub_CONUS <- readRDS("./models/betaLASSO/ShrubCover_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
# Broad-leaved tree:
# Grass/shrub: best lambda model
BLtree_GS <- readRDS("./models/betaLASSO/AngioTreeCover_prop_shrubGrass_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
# Forest: best lambda model
BLtree_F <- readRDS("./models/betaLASSO/AngioTreeCover_prop_forest_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
# Needle-leaved tree:
# Grass/shrub: best lambda model
NLtree_GS <- readRDS("./models/betaLASSO/ConifTreeCover_prop_shrubGrass_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
# Forest: best lambda model
NLtree_F <- readRDS("./models/betaLASSO/ConifTreeCover_prop_forest_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
# Forb:
# CONUS: best lambda model
forb_CONUS <- readRDS("./models/betaLASSO/ForbCover_prop_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
# C3 grass:
# CONUS: 1SE lambda model
C3grass_CONUS <- readRDS("./models/betaLASSO/C3GramCover_prop_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_oneSELambdaGLM.rds")
# C4 grass:
# CONUS: best lambda model
C4grass_CONUS <- readRDS("./models/betaLASSO/C4GramCover_prop_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")Now, predict absolute cover with scaled climate/weather/soils data
## with contemporary climate data
totalHerb_F_predsContemp <- makePredictions(predictionDF = climDatPred,
modelObject = totalHerb_F) %>%
rename(absTotalHerb_F = "modelPreds")
totalHerb_GS_predsContemp <- makePredictions(predictionDF = climDatPred,
modelObject = totalHerb_GS) %>%
rename(absTotalHerb_GS = "modelPreds")
totalTree_F_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = totalTree_F) %>%
rename(absTotalTree_F = "modelPreds")
totalTree_GS_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = totalTree_GS) %>%
rename(absTotalTree_GS = "modelPreds")
bareGround_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = bareGround_CONUS) %>%
rename(absBareGround_CONUS = "modelPreds")
shrub_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = shrub_CONUS) %>%
rename(absShrub_CONUS = "modelPreds")
BLtree_F_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = BLtree_F) %>%
rename(percBLTree_F = "modelPreds")
BLtree_GS_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = BLtree_GS) %>%
rename(percBLTree_GS = "modelPreds")
NLtree_F_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = NLtree_F) %>%
rename(percNLTree_F = "modelPreds")
NLtree_GS_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = NLtree_GS) %>%
rename(percNLTree_GS = "modelPreds")
C3grass_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = C3grass_CONUS) %>%
rename(percC3grass_CONUS = "modelPreds")
C4grass_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = C4grass_CONUS) %>%
rename(percC4grass_CONUS = "modelPreds")
forb_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = forb_CONUS) %>%
rename(percForb_CONUS = "modelPreds")
contempPreds <- totalHerb_F_predsContemp %>%
cbind(totalHerb_GS_predsContemp %>% select(absTotalHerb_GS)) %>%
cbind(totalTree_F_predsContemp %>% select(absTotalTree_F)) %>%
cbind(totalTree_GS_predsContemp %>% select(absTotalTree_GS)) %>%
cbind(bareGround_CONUS_predsContemp %>% select(absBareGround_CONUS)) %>%
cbind(shrub_CONUS_predsContemp %>% select(absShrub_CONUS)) %>%
cbind(BLtree_F_predsContemp %>% select(percBLTree_F)) %>%
cbind(BLtree_GS_predsContemp %>% select(percBLTree_GS)) %>%
cbind(NLtree_F_predsContemp %>% select(percNLTree_F)) %>%
cbind(NLtree_GS_predsContemp %>% select(percNLTree_GS)) %>%
cbind(C3grass_CONUS_predsContemp %>% select(percC3grass_CONUS)) %>%
cbind(C4grass_CONUS_predsContemp %>% select(percC4grass_CONUS)) %>%
cbind(forb_CONUS_predsContemp %>% select(percForb_CONUS)) %>%
cbind(preds_byHand %>% select(pred)) %>%
rename(prob_Forest = "pred") %>%
mutate(prob_grassShrub = 1 - prob_Forest)
# get predictions with first climate model
## with contemporary climate data
totalHerb_F_predsFuture1 <- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = totalHerb_F) %>%
rename(absTotalHerb_F = "modelPreds")
totalHerb_GS_predsFuture1 <- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = totalHerb_GS) %>%
rename(absTotalHerb_GS = "modelPreds")
totalTree_F_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = totalTree_F) %>%
rename(absTotalTree_F = "modelPreds")
totalTree_GS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = totalTree_GS) %>%
rename(absTotalTree_GS = "modelPreds")
bareGround_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = bareGround_CONUS) %>%
rename(absBareGround_CONUS = "modelPreds")
shrub_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = shrub_CONUS) %>%
rename(absShrub_CONUS = "modelPreds")
BLtree_F_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = BLtree_F) %>%
rename(percBLTree_F = "modelPreds")
BLtree_GS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = BLtree_GS) %>%
rename(percBLTree_GS = "modelPreds")
NLtree_F_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = NLtree_F) %>%
rename(percNLTree_F = "modelPreds")
NLtree_GS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = NLtree_GS) %>%
rename(percNLTree_GS = "modelPreds")
C3grass_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = C3grass_CONUS) %>%
rename(percC3grass_CONUS = "modelPreds")
C4grass_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = C4grass_CONUS) %>%
rename(percC4grass_CONUS = "modelPreds")
forb_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = forb_CONUS) %>%
rename(percForb_CONUS = "modelPreds")
predsFuture1 <- totalHerb_F_predsFuture1 %>%
cbind(totalHerb_GS_predsFuture1 %>% select(absTotalHerb_GS)) %>%
cbind(totalTree_F_predsFuture1 %>% select(absTotalTree_F)) %>%
cbind(totalTree_GS_predsFuture1 %>% select(absTotalTree_GS)) %>%
cbind(bareGround_CONUS_predsFuture1 %>% select(absBareGround_CONUS)) %>%
cbind(shrub_CONUS_predsFuture1 %>% select(absShrub_CONUS)) %>%
cbind(BLtree_F_predsFuture1 %>% select(percBLTree_F)) %>%
cbind(BLtree_GS_predsFuture1 %>% select(percBLTree_GS)) %>%
cbind(NLtree_F_predsFuture1 %>% select(percNLTree_F)) %>%
cbind(NLtree_GS_predsFuture1 %>% select(percNLTree_GS)) %>%
cbind(C3grass_CONUS_predsFuture1 %>% select(percC3grass_CONUS)) %>%
cbind(C4grass_CONUS_predsFuture1 %>% select(percC4grass_CONUS)) %>%
cbind(forb_CONUS_predsFuture1 %>% select(percForb_CONUS)) %>%
cbind(preds_future1_byHand %>% select(pred)) %>%
rename(prob_Forest = "pred") %>%
mutate(prob_grassShrub = 1 - prob_Forest)
# get predictions with second climate model
totalHerb_F_predsFuture2 <- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = totalHerb_F) %>%
rename(absTotalHerb_F = "modelPreds")
totalHerb_GS_predsFuture2 <- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = totalHerb_GS) %>%
rename(absTotalHerb_GS = "modelPreds")
totalTree_F_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = totalTree_F) %>%
rename(absTotalTree_F = "modelPreds")
totalTree_GS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = totalTree_GS) %>%
rename(absTotalTree_GS = "modelPreds")
bareGround_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = bareGround_CONUS) %>%
rename(absBareGround_CONUS = "modelPreds")
shrub_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = shrub_CONUS) %>%
rename(absShrub_CONUS = "modelPreds")
BLtree_F_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = BLtree_F) %>%
rename(percBLTree_F = "modelPreds")
BLtree_GS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = BLtree_GS) %>%
rename(percBLTree_GS = "modelPreds")
NLtree_F_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = NLtree_F) %>%
rename(percNLTree_F = "modelPreds")
NLtree_GS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = NLtree_GS) %>%
rename(percNLTree_GS = "modelPreds")
C3grass_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = C3grass_CONUS) %>%
rename(percC3grass_CONUS = "modelPreds")
C4grass_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = C4grass_CONUS) %>%
rename(percC4grass_CONUS = "modelPreds")
forb_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = forb_CONUS) %>%
rename(percForb_CONUS = "modelPreds")
predsFuture2 <- totalHerb_F_predsFuture2 %>%
cbind(totalHerb_GS_predsFuture2 %>% select(absTotalHerb_GS)) %>%
cbind(totalTree_F_predsFuture2 %>% select(absTotalTree_F)) %>%
cbind(totalTree_GS_predsFuture2 %>% select(absTotalTree_GS)) %>%
cbind(bareGround_CONUS_predsFuture2 %>% select(absBareGround_CONUS)) %>%
cbind(shrub_CONUS_predsFuture2 %>% select(absShrub_CONUS)) %>%
cbind(BLtree_F_predsFuture2 %>% select(percBLTree_F)) %>%
cbind(BLtree_GS_predsFuture2 %>% select(percBLTree_GS)) %>%
cbind(NLtree_F_predsFuture2 %>% select(percNLTree_F)) %>%
cbind(NLtree_GS_predsFuture2 %>% select(percNLTree_GS)) %>%
cbind(C3grass_CONUS_predsFuture2 %>% select(percC3grass_CONUS)) %>%
cbind(C4grass_CONUS_predsFuture2 %>% select(percC4grass_CONUS)) %>%
cbind(forb_CONUS_predsFuture2 %>% select(percForb_CONUS)) %>%
cbind(preds_future2_byHand %>% select(pred)) %>%
rename(prob_Forest = "pred") %>%
mutate(prob_grassShrub = 1 - prob_Forest)“Scale” absolute cover predictions according to predicted ecoregion probability
## for contemporary data
contempPreds <- contempPreds %>%
mutate(absTotalHerb_CONUS = (prob_grassShrub * absTotalHerb_GS) + (prob_Forest * absTotalHerb_F),
absTotalTree_CONUS = (prob_grassShrub * absTotalTree_GS) + (prob_Forest * absTotalTree_F),
percBLTree_CONUS = (prob_grassShrub * percBLTree_GS) + (prob_Forest * percBLTree_F),
percNLTree_CONUS = (prob_grassShrub * percNLTree_GS) + (prob_Forest * percNLTree_F)
)
## for first climate model
predsFuture1 <- predsFuture1 %>%
mutate(absTotalHerb_CONUS = (prob_grassShrub * absTotalHerb_GS) + (prob_Forest * absTotalHerb_F),
absTotalTree_CONUS = (prob_grassShrub * absTotalTree_GS) + (prob_Forest * absTotalTree_F),
percBLTree_CONUS = (prob_grassShrub * percBLTree_GS) + (prob_Forest * percBLTree_F),
percNLTree_CONUS = (prob_grassShrub * percNLTree_GS) + (prob_Forest * percNLTree_F)
)
## for second climate model
predsFuture2 <- predsFuture2 %>%
mutate(absTotalHerb_CONUS = (prob_grassShrub * absTotalHerb_GS) + (prob_Forest * absTotalHerb_F),
absTotalTree_CONUS = (prob_grassShrub * absTotalTree_GS) + (prob_Forest * absTotalTree_F),
percBLTree_CONUS = (prob_grassShrub * percBLTree_GS) + (prob_Forest * percBLTree_F),
percNLTree_CONUS = (prob_grassShrub * percNLTree_GS) + (prob_Forest * percNLTree_F)
)For percentages (NL/BL tree and C4/C3/Forb), scale so that they sum to one
contempPreds <- contempPreds %>%
mutate(percBLTree_scaled_CONUS = percBLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
percNLTree_scaled_CONUS = percNLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
percC3grass_scaled_CONUS = percC3grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
percC4grass_scaled_CONUS = percC4grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
percForb_scaled_CONUS = percForb_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS))
predsFuture1 <- predsFuture1 %>%
mutate(percBLTree_scaled_CONUS = percBLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
percNLTree_scaled_CONUS = percNLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
percC3grass_scaled_CONUS = percC3grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
percC4grass_scaled_CONUS = percC4grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
percForb_scaled_CONUS = percForb_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS))
predsFuture2 <- predsFuture2 %>%
mutate(percBLTree_scaled_CONUS = percBLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
percNLTree_scaled_CONUS = percNLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
percC3grass_scaled_CONUS = percC3grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
percC4grass_scaled_CONUS = percC4grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
percForb_scaled_CONUS = percForb_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS))Convert percentages for level 2 cover into absolute cover values
## for contemporary data
contempPreds <- contempPreds %>%
mutate(absNLTree_CONUS = (percNLTree_scaled_CONUS * absTotalTree_CONUS),
absBLTree_CONUS = (percBLTree_scaled_CONUS * absTotalTree_CONUS),
absC3grass_CONUS = (percC3grass_scaled_CONUS * absTotalHerb_CONUS),
absC4grass_CONUS = (percC4grass_scaled_CONUS * absTotalHerb_CONUS),
absForb_CONUS = (percForb_scaled_CONUS * absTotalHerb_CONUS)
)
## for first climate model
predsFuture1 <- predsFuture1 %>%
mutate(absNLTree_CONUS = (percNLTree_scaled_CONUS * absTotalTree_CONUS),
absBLTree_CONUS = (percBLTree_scaled_CONUS * absTotalTree_CONUS),
absC3grass_CONUS = (percC3grass_scaled_CONUS * absTotalHerb_CONUS),
absC4grass_CONUS = (percC4grass_scaled_CONUS * absTotalHerb_CONUS),
absForb_CONUS = (percForb_scaled_CONUS * absTotalHerb_CONUS)
)
## for second climate model
predsFuture2 <- predsFuture2 %>%
mutate(absNLTree_CONUS = (percNLTree_scaled_CONUS * absTotalTree_CONUS),
absBLTree_CONUS = (percBLTree_scaled_CONUS * absTotalTree_CONUS),
absC3grass_CONUS = (percC3grass_scaled_CONUS * absTotalHerb_CONUS),
absC4grass_CONUS = (percC4grass_scaled_CONUS * absTotalHerb_CONUS),
absForb_CONUS = (percForb_scaled_CONUS * absTotalHerb_CONUS)
)Show predictions of absolute cover, CONUS-wide (after scaling according to ecoregion classification probability)
# prep observed data for use in figures
obsDat <- modDat_1_s %>%
select(Long, Lat, Year, ShrubCover, TotalHerbaceousCover, TotalTreeCover, BareGroundCover, C3GramCover_prop, C4GramCover_prop, ForbCover_prop, AngioTreeCover_prop, ConifTreeCover_prop) %>%
mutate(absNLTree_CONUS = ConifTreeCover_prop * TotalTreeCover,
absBLTree_CONUS = AngioTreeCover_prop * TotalTreeCover,
absC3grass_CONUS = C3GramCover_prop * TotalTreeCover,
absC4grass_CONUS = C4GramCover_prop * TotalTreeCover,
absForb_CONUS = ForbCover_prop * TotalTreeCover
) %>%
rename(absShrub_CONUS = ShrubCover,
absTotalTree_CONUS = TotalTreeCover,
absTotalHerb_CONUS = TotalHerbaceousCover,
absBareGround_CONUS = BareGroundCover)
# list of absolute cover names
responseNames <- c("absTotalTree_CONUS", "absTotalHerb_CONUS", "absShrub_CONUS", "absBareGround_CONUS", "absNLTree_CONUS", "absBLTree_CONUS", "absC3grass_CONUS", "absC4grass_CONUS", "absForb_CONUS")
absoluteCoverMaps_contemp <- lapply(responseNames,
FUN = function(x) {
## contemporary absolute cover
# rasterize response
resp_rast <- contempPreds %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = x,
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
resp_rast_2 <- resp_rast %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# make a map
resp_map <- ggplot() +
geom_spatraster(data = resp_rast_2) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Predictions of absolute cover of ",x,",
\n after scaling according to ecoregion prediction")) +
scale_fill_gradient2(low = "brown",
mid = "wheat" ,
high = "darkgreen" ,
midpoint = 0, limits = c(0,1), na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist <- ggplot() +
geom_density(aes(contempPreds[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of absolute cover of ",x)) +
theme_classic()
## calculate residuals
# make maps of residuals between absolute cover predictions and observations
# rasterize observations
plotObservations <- obsDat %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = x,
fun = mean, na.rm = TRUE)
plotObservations <- plotObservations %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# calculate residuals (observed - predicted)
resids_rast <- plotObservations - resp_rast_2
## aggregate the residuals raster to make it easier to see
resids_rastCoarse <- resids_rast %>%
aggregate(fact = 2, fun = "mean")
# map the residuals
mapResids <- ggplot() +
geom_spatraster(data = resids_rastCoarse) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(resids_rast)),fill=NA ) +
labs(title = paste0("Residuals (observations - predicted absolute cover) for ",x)) +
scale_fill_gradient2(low = "red",
mid = "white" ,
high = "blue" ,
midpoint = 0, na.value = "lightgrey",
limits = c(-1,1)
) +
xlim(st_bbox(resids_rast)[c(1,3)]) +
ylim(st_bbox(resids_rast)[c(2,4)]) +
theme_classic()
# make a histogram of residuals
histResids <- ggplot() +
geom_density(aes(terra::values(resids_rast$mean, na.rm = TRUE, mat = FALSE)),
col = "black", fill = "darkgrey") +
geom_vline(aes(xintercept = 0), linetype = 2) +
xlab(paste0("Residuals (obs - pred) ",x)) +
theme_classic()
## future model 1 absolute cover
# rasterize response
resp_rast_future1 <- predsFuture1 %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = x,
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
resp_rast_future1_2 <- resp_rast_future1 %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# calculate differences between contemporary preds and future 1 preds (deltas)
rast_deltas_future1 <- resp_rast_future1_2 - resp_rast_2
# make a map
resp_map_future1 <- ggplot() +
geom_spatraster(data = rast_deltas_future1) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Future predictions of absolute cover of ",x,", \n compared to contemporary predictions (future - contemprary) \n - climate model: BNU-ESM")) +
scale_fill_gradient2(low = "orange",
mid = "white" ,
high = "purple" ,
midpoint = 0, limits = c(-1,1), na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist_future1 <- ggplot() +
geom_density(aes(predsFuture1[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of absolute cover of ",x,
" \n in future scenario (model = BNU-ESM)")) +
theme_classic()
## future model 2 absolute cover
# rasterize response
resp_rast_future2 <- predsFuture2 %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = x,
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
resp_rast_future2_2 <- resp_rast_future2 %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# calculate differences between contemporary preds and future 1 preds (deltas)
rast_deltas_future2 <- resp_rast_future2_2 - resp_rast_2
# make a map
resp_map_future2 <- ggplot() +
geom_spatraster(data = rast_deltas_future2) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Future predictions of absolute cover of ",x,", \n compared to contemporary predictions (future - contemprary) \n- climate model: IPSL-CM5A-MR (France)")) +
scale_fill_gradient2(low = "orange",
mid = "white" ,
high = "purple" ,
midpoint = 0, limits = c(-1,1), na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist_future2 <- ggplot() +
geom_density(aes(predsFuture2[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of absolute cover of ",x,
" \n in future scenario (model = IPSL-CM5A-MR (France))")) +
theme_classic()
return(list("rast" = resp_rast_2,
"map" = resp_map,
"hist" = resp_hist,
"rast_resids" = resids_rast,
"map_resids" = mapResids,
"hist_resids" = histResids,
"rast_future1" = resp_rast_future1_2,
"map_future1" = resp_map_future1,
"hist_future1" = resp_hist_future1,
"rast_future2" = resp_rast_future2_2,
"map_future2" = resp_map_future2,
"hist_future2" = resp_hist_future2))
})
names(absoluteCoverMaps_contemp) <- c("absTotalTree_CONUS", "absTotalHerb_CONUS", "absShrub_CONUS", "absBareGround_CONUS", "absNLTree_CONUS", "absBLTree_CONUS", "absC3grass_CONUS", "absC4grass_CONUS", "absForb_CONUS")ggarrange(
ggarrange(absoluteCoverMaps_contemp$absTotalTree_CONUS$map,
absoluteCoverMaps_contemp$absTotalTree_CONUS$map_resids,
absoluteCoverMaps_contemp$absTotalTree_CONUS$map_future1,
absoluteCoverMaps_contemp$absTotalTree_CONUS$map_future2,
absoluteCoverMaps_contemp$absTotalTree_CONUS$hist,
absoluteCoverMaps_contemp$absTotalTree_CONUS$hist_resids,
absoluteCoverMaps_contemp$absTotalTree_CONUS$hist_future1,
absoluteCoverMaps_contemp$absTotalTree_CONUS$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(absoluteCoverMaps_contemp$absTotalHerb_CONUS$map,
absoluteCoverMaps_contemp$absTotalHerb_CONUS$map_resids,
absoluteCoverMaps_contemp$absTotalHerb_CONUS$map_future1,
absoluteCoverMaps_contemp$absTotalHerb_CONUS$map_future2,
absoluteCoverMaps_contemp$absTotalHerb_CONUS$hist,
absoluteCoverMaps_contemp$absTotalHerb_CONUS$hist_resids,
absoluteCoverMaps_contemp$absTotalHerb_CONUS$hist_future1,
absoluteCoverMaps_contemp$absTotalHerb_CONUS$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(absoluteCoverMaps_contemp$absShrub_CONUS$map,
absoluteCoverMaps_contemp$absShrub_CONUS$map_resids,
absoluteCoverMaps_contemp$absShrub_CONUS$map_future1,
absoluteCoverMaps_contemp$absShrub_CONUS$map_future2,
absoluteCoverMaps_contemp$absShrub_CONUS$hist,
absoluteCoverMaps_contemp$absShrub_CONUS$hist_resids,
absoluteCoverMaps_contemp$absShrub_CONUS$hist_future1,
absoluteCoverMaps_contemp$absShrub_CONUS$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(absoluteCoverMaps_contemp$absBareGround_CONUS$map,
absoluteCoverMaps_contemp$absBareGround_CONUS$map_resids,
absoluteCoverMaps_contemp$absBareGround_CONUS$map_future1,
absoluteCoverMaps_contemp$absBareGround_CONUS$map_future2,
absoluteCoverMaps_contemp$absBareGround_CONUS$hist,
absoluteCoverMaps_contemp$absBareGround_CONUS$hist_resids,
absoluteCoverMaps_contemp$absBareGround_CONUS$hist_future1,
absoluteCoverMaps_contemp$absBareGround_CONUS$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(absoluteCoverMaps_contemp$absNLTree_CONUS$map,
absoluteCoverMaps_contemp$absNLTree_CONUS$map_resids,
absoluteCoverMaps_contemp$absNLTree_CONUS$map_future1,
absoluteCoverMaps_contemp$absNLTree_CONUS$map_future2,
absoluteCoverMaps_contemp$absNLTree_CONUS$hist,
absoluteCoverMaps_contemp$absNLTree_CONUS$hist_resids,
absoluteCoverMaps_contemp$absNLTree_CONUS$hist_future1,
absoluteCoverMaps_contemp$absNLTree_CONUS$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(absoluteCoverMaps_contemp$absBLTree_CONUS$map,
absoluteCoverMaps_contemp$absBLTree_CONUS$map_resids,
absoluteCoverMaps_contemp$absBLTree_CONUS$map_future1,
absoluteCoverMaps_contemp$absBLTree_CONUS$map_future2,
absoluteCoverMaps_contemp$absBLTree_CONUS$hist,
absoluteCoverMaps_contemp$absBLTree_CONUS$hist_resids,
absoluteCoverMaps_contemp$absBLTree_CONUS$hist_future1,
absoluteCoverMaps_contemp$absBLTree_CONUS$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(absoluteCoverMaps_contemp$absC3grass_CONUS$map,
absoluteCoverMaps_contemp$absC3grass_CONUS$map_resids,
absoluteCoverMaps_contemp$absC3grass_CONUS$map_future1,
absoluteCoverMaps_contemp$absC3grass_CONUS$map_future2,
absoluteCoverMaps_contemp$absC3grass_CONUS$hist,
absoluteCoverMaps_contemp$absC3grass_CONUS$hist_resids,
absoluteCoverMaps_contemp$absC3grass_CONUS$hist_future1,
absoluteCoverMaps_contemp$absC3grass_CONUS$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(absoluteCoverMaps_contemp$absC4grass_CONUS$map,
absoluteCoverMaps_contemp$absC4grass_CONUS$map_resids,
absoluteCoverMaps_contemp$absC4grass_CONUS$map_future1,
absoluteCoverMaps_contemp$absC4grass_CONUS$map_future2,
absoluteCoverMaps_contemp$absC4grass_CONUS$hist,
absoluteCoverMaps_contemp$absC4grass_CONUS$hist_resids,
absoluteCoverMaps_contemp$absC4grass_CONUS$hist_future1,
absoluteCoverMaps_contemp$absC4grass_CONUS$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(absoluteCoverMaps_contemp$absForb_CONUS$map,
absoluteCoverMaps_contemp$absForb_CONUS$map_resids,
absoluteCoverMaps_contemp$absForb_CONUS$map_future1,
absoluteCoverMaps_contemp$absForb_CONUS$map_future2,
absoluteCoverMaps_contemp$absForb_CONUS$hist,
absoluteCoverMaps_contemp$absForb_CONUS$hist_resids,
absoluteCoverMaps_contemp$absForb_CONUS$hist_future1,
absoluteCoverMaps_contemp$absForb_CONUS$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ncol = 1
)
Get predictions of absolute cover across time using observed climate/weather/soils data
# get climate data from dayMet (d.f. is "climDatPred_unscaled")
climDat_time_unscaled <- climDat_time
names(climDat_time_unscaled)[c(5, 6, 7, 13, 10, 12, 99, 101, 102, 103)] <- c("T_warmestMonth_meanAnnAvg_CLIM", "T_coldestMonth_meanAnnAvg_CLIM", "precip_wettestMonth_meanAnnAvg_CLIM", "annWaterDeficit_meanAnnAvg_CLIM", "PrecipTempCorr_meanAnnAvg_CLIM", "isothermality_meanAnnAvg_CLIM", "soilDepth", "avgSandPerc_acrossDepth", "avgCoarsePerc_acrossDepth", "avgOrganicCarbonPerc_0_3cm")
# predict with contemporary climate data
preds_byHand_time <- climDat_time_unscaled %>%
mutate(pred = 1/(1 + exp(-( 9.8726 + -0.2999*T_warmestMonth_meanAnnAvg_CLIM + 0.2456*T_coldestMonth_meanAnnAvg_CLIM + 0.0106*precip_wettestMonth_meanAnnAvg_CLIM + -0.0621*annWaterDeficit_meanAnnAvg_CLIM + -2.7863*PrecipTempCorr_meanAnnAvg_CLIM + 0.0540*isothermality_meanAnnAvg_CLIM + -0.0076*soilDepth + 0.0335*avgSandPerc_acrossDepth + 0.0310*avgCoarsePerc_acrossDepth + 0.2726*avgOrganicCarbonPerc_0_3cm))))
#predictionsModel <- predict(object = ecoMod, type = "response", newdata = climDatPred_unscaled)
## with contemporary climate data
totalHerb_F_predsOverTime <- makePredictions(predictionDF = climDat_time_Pred,
modelObject = totalHerb_F) %>%
rename(absTotalHerb_F = "modelPreds")
totalHerb_GS_predsOverTime <- makePredictions(predictionDF = climDat_time_Pred,
modelObject = totalHerb_GS) %>%
rename(absTotalHerb_GS = "modelPreds")
totalTree_F_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = totalTree_F) %>%
rename(absTotalTree_F = "modelPreds")
totalTree_GS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = totalTree_GS) %>%
rename(absTotalTree_GS = "modelPreds")
bareGround_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = bareGround_CONUS) %>%
rename(absBareGround_CONUS = "modelPreds")
shrub_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = shrub_CONUS) %>%
rename(absShrub_CONUS = "modelPreds")
BLtree_F_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = BLtree_F) %>%
rename(percBLTree_F = "modelPreds")
BLtree_GS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = BLtree_GS) %>%
rename(percBLTree_GS = "modelPreds")
NLtree_F_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = NLtree_F) %>%
rename(percNLTree_F = "modelPreds")
NLtree_GS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = NLtree_GS) %>%
rename(percNLTree_GS = "modelPreds")
C3grass_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = C3grass_CONUS) %>%
rename(percC3grass_CONUS = "modelPreds")
C4grass_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = C4grass_CONUS) %>%
rename(percC4grass_CONUS = "modelPreds")
forb_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = forb_CONUS) %>%
rename(percForb_CONUS = "modelPreds")
predsOverTime <- totalHerb_F_predsOverTime %>%
cbind(totalHerb_GS_predsOverTime %>% select(absTotalHerb_GS)) %>%
cbind(totalTree_F_predsOverTime %>% select(absTotalTree_F)) %>%
cbind(totalTree_GS_predsOverTime %>% select(absTotalTree_GS)) %>%
cbind(bareGround_CONUS_predsOverTime %>% select(absBareGround_CONUS)) %>%
cbind(shrub_CONUS_predsOverTime %>% select(absShrub_CONUS)) %>%
cbind(BLtree_F_predsOverTime %>% select(percBLTree_F)) %>%
cbind(BLtree_GS_predsOverTime %>% select(percBLTree_GS)) %>%
cbind(NLtree_F_predsOverTime %>% select(percNLTree_F)) %>%
cbind(NLtree_GS_predsOverTime %>% select(percNLTree_GS)) %>%
cbind(C3grass_CONUS_predsOverTime %>% select(percC3grass_CONUS)) %>%
cbind(C4grass_CONUS_predsOverTime %>% select(percC4grass_CONUS)) %>%
cbind(forb_CONUS_predsOverTime %>% select(percForb_CONUS)) %>%
cbind(preds_byHand_time %>% select(pred)) %>%
rename(prob_Forest = "pred") %>%
mutate(prob_grassShrub = 1 - prob_Forest)
## now, scale predictions according to ecoregion percentage
predsOverTime <- predsOverTime %>%
mutate(absTotalHerb_CONUS = (prob_grassShrub * absTotalHerb_GS) + (prob_Forest * absTotalHerb_F),
absTotalTree_CONUS = (prob_grassShrub * absTotalTree_GS) + (prob_Forest * absTotalTree_F),
percBLTree_CONUS = (prob_grassShrub * percBLTree_GS) + (prob_Forest * percBLTree_F),
percNLTree_CONUS = (prob_grassShrub * percNLTree_GS) + (prob_Forest * percNLTree_F)
)
## then, scale the percentages so they sum to 1
predsOverTime <- predsOverTime %>%
mutate(percBLTree_scaled_CONUS = percBLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
percNLTree_scaled_CONUS = percNLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
percC3grass_scaled_CONUS = percC3grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
percC4grass_scaled_CONUS = percC4grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
percForb_scaled_CONUS = percForb_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS))
## convert percentages into absolute cover for level 2 cover classes
predsOverTime <- predsOverTime %>%
mutate(absNLTree_CONUS = (percNLTree_scaled_CONUS * absTotalTree_CONUS),
absBLTree_CONUS = (percBLTree_scaled_CONUS * absTotalTree_CONUS),
absC3grass_CONUS = (percC3grass_scaled_CONUS * absTotalHerb_CONUS),
absC4grass_CONUS = (percC4grass_scaled_CONUS * absTotalHerb_CONUS),
absForb_CONUS = (percForb_scaled_CONUS * absTotalHerb_CONUS)
)Now, visualize predictions of absolute cover over time
# change spatialIDs to discrete numbers, since the long form is causing some sort of floating point issue?
spatID_lu <- data.frame("spatialID" = unique(predsOverTime$spatialID),
"uniqueID" = 1:length(predsOverTime$spatialID))
predsOverTime <- predsOverTime %>%
left_join(spatID_lu)
# make into a long format
predsOverTime_long <- predsOverTime %>%
select(year,x, y, spatialID, uniqueID, absShrub_CONUS, absBareGround_CONUS, absTotalHerb_CONUS, absTotalTree_CONUS, absNLTree_CONUS:absForb_CONUS) %>%
pivot_longer(cols = c(absShrub_CONUS, absBareGround_CONUS, absTotalHerb_CONUS, absTotalTree_CONUS, absNLTree_CONUS:absForb_CONUS)) %>%
mutate(year = as.integer(year))
# calculate the coefficient of variation for each cover type at each location
predsOverTime_cv <- predsOverTime_long %>%
group_by(uniqueID, name) %>%
summarize(cv = sd(value)/mean(value))
## add x,y
predsOverTime_cv <- predsOverTime_cv %>%
left_join(predsOverTime %>% select(uniqueID, x, y))
# visualize change over time for a selection of points
predsOverTime_long %>%
filter(uniqueID %in% 1:50) %>%
filter(name %in% c("absShrub_CONUS", "absBareGround_CONUS", "absNLTree_CONUS",
"absBLTree_CONUS", "absC3grass_CONUS", "absC4grass_CONUS", "absForb_CONUS")) %>%
ggplot() +
facet_wrap(~uniqueID) +
geom_area(aes(x = year, y = value, fill = name)) +
geom_line(data = predsOverTime_long %>%
filter(uniqueID %in% 1:50) %>%
filter(name %in% c("absTotalHerb_CONUS", "absTotalTree_CONUS")),
aes(x = year, y = value, linetype = name), col = "black") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90)) +
labs(title = "Change in modeled absolute cover by functional type from \n 2010 to 2020 for a subset of locations")
predsOverTime_cv %>%
ggplot() +
geom_boxplot(aes(x = name, y = cv, col = name)) +
theme_minimal() +
labs(title = "Coefficients of variation of model-predicted absolute \n cover at a given location for each functional group from 2010-2020",
x = "functional group",
y = "coefficient of variation") +
theme(axis.text.x = element_text(angle = 90))
# show the locations of each point that is being used for this temporal analysis
#cropped_states <- cropped_states %>%
# rename(geometry = Shape)
ggplot() +
facet_wrap(~name) +
geom_sf(data = cropped_states, aes()) +
stat_summary2d(data = predsOverTime_cv, aes(x = x, y = y, z = cv), fun = ~ mean(.x, na.rm = TRUE)) +
geom_hex(data = predsOverTime_cv, aes(x = x, y = y, col = cv), stat = "mean", na.rm = TRUE) +
guides(fill = guide_legend(title = "CV")) +
theme_minimal() +
labs(title = "Locations of points included in the random \n sample for examining temporal variation \n color-coded by the coefficient of variation ")
Now, relativize the predictions across the different functional groups
First option: relativize cover so that cover for all groups sums to 1
## for contemporary data
contempPreds <- contempPreds %>%
mutate(totalAbsCover_CONUS = (absTotalHerb_CONUS + absTotalTree_CONUS + absBareGround_CONUS + absShrub_CONUS)) %>%
mutate(relCoverA_totalTree = absTotalTree_CONUS/totalAbsCover_CONUS,
relCoverA_totalHerb = absTotalHerb_CONUS/totalAbsCover_CONUS,
relCoverA_shrub = absShrub_CONUS/totalAbsCover_CONUS,
relCoverA_bareGround = absBareGround_CONUS/totalAbsCover_CONUS,
relCoverA_NLTree = absNLTree_CONUS/totalAbsCover_CONUS,
relCoverA_BLTree = absBLTree_CONUS/totalAbsCover_CONUS,
relCoverA_C3Grass = absC3grass_CONUS/totalAbsCover_CONUS,
relCoverA_C4Grass = absC4grass_CONUS/totalAbsCover_CONUS,
relCoverA_Forb = absForb_CONUS/totalAbsCover_CONUS
)
## for first climate model
predsFuture1 <- predsFuture1 %>%
mutate(totalAbsCover_CONUS = (absTotalHerb_CONUS + absTotalTree_CONUS + absBareGround_CONUS + absShrub_CONUS)) %>%
mutate(relCoverA_totalTree = absTotalTree_CONUS/totalAbsCover_CONUS,
relCoverA_totalHerb = absTotalHerb_CONUS/totalAbsCover_CONUS,
relCoverA_shrub = absShrub_CONUS/totalAbsCover_CONUS,
relCoverA_bareGround = absBareGround_CONUS/totalAbsCover_CONUS,
relCoverA_NLTree = absNLTree_CONUS/totalAbsCover_CONUS,
relCoverA_BLTree = absBLTree_CONUS/totalAbsCover_CONUS,
relCoverA_C3Grass = absC3grass_CONUS/totalAbsCover_CONUS,
relCoverA_C4Grass = absC4grass_CONUS/totalAbsCover_CONUS,
relCoverA_Forb = absForb_CONUS/totalAbsCover_CONUS
)
## for second climate model
predsFuture2 <- predsFuture2 %>%
mutate(totalAbsCover_CONUS = (absTotalHerb_CONUS + absTotalTree_CONUS + absBareGround_CONUS + absShrub_CONUS)) %>%
mutate(relCoverA_totalTree = absTotalTree_CONUS/totalAbsCover_CONUS,
relCoverA_totalHerb = absTotalHerb_CONUS/totalAbsCover_CONUS,
relCoverA_shrub = absShrub_CONUS/totalAbsCover_CONUS,
relCoverA_bareGround = absBareGround_CONUS/totalAbsCover_CONUS,
relCoverA_NLTree = absNLTree_CONUS/totalAbsCover_CONUS,
relCoverA_BLTree = absBLTree_CONUS/totalAbsCover_CONUS,
relCoverA_C3Grass = absC3grass_CONUS/totalAbsCover_CONUS,
relCoverA_C4Grass = absC4grass_CONUS/totalAbsCover_CONUS,
relCoverA_Forb = absForb_CONUS/totalAbsCover_CONUS
)Now, show the results in map form
# list of absolute cover names
responseNames <- c("relCoverA_totalTree", "relCoverA_totalHerb", "relCoverA_shrub", "relCoverA_bareGround",
"relCoverA_NLTree", "relCoverA_BLTree", "relCoverA_C3Grass", "relCoverA_C4Grass", "relCoverA_Forb",
"totalAbsCover_CONUS")
relCover_A_Maps_contemp <- lapply(responseNames,
FUN = function(x) {
## contemporary absolute cover
# rasterize response
resp_rast <- contempPreds %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = x,
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
resp_rast_2 <- resp_rast %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
## future model 1 absolute cover
# rasterize response
resp_rast_future1 <- predsFuture1 %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = x,
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
resp_rast_future1_2 <- resp_rast_future1 %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# calculate differences between contemporary preds and future 1 preds (deltas)
rast_deltas_future1 <- resp_rast_future1_2 - resp_rast_2
## future model 2 absolute cover
# rasterize response
resp_rast_future2 <- predsFuture2 %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = x,
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
resp_rast_future2_2 <- resp_rast_future2 %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# calculate differences between contemporary preds and future 1 preds (deltas)
rast_deltas_future2 <- resp_rast_future2_2 - resp_rast_2
if (x == "totalAbsCover_CONUS") {
# make a map of predictions w/ current data
resp_map <- ggplot() +
geom_spatraster(data = resp_rast_2) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Predictions of Total Absolute Cover with contemporary data")) +
scale_fill_gradient2(low = "lightblue",
#mid = "wheat" ,
high = "darkblue" ,
#midpoint = 0,
limits = c(0,max(terra::values(resp_rast_2$mean), na.rm= TRUE)),
na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist <- ggplot() +
geom_density(aes(contempPreds[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of Total Absolute Cover with contemporary data")) +
theme_classic()
# make a map of future preds w/ model 1
resp_map_future1 <- ggplot() +
geom_spatraster(data = rast_deltas_future1) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Future predictions of Total Absolute Cover \n compared to contemporary predictions (future - contemprary) \n - climate model: BNU-ESM")) +
scale_fill_gradient2(low = "orange",
mid = "white" ,
high = "purple" ,
midpoint = 0, limits = c(-2,2), na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist_future1 <- ggplot() +
geom_density(aes(predsFuture1[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of relative cover of ",x,
" \n in future scenario (model = BNU-ESM)")) +
theme_classic()
# make a map of future preds w/ model 2
resp_map_future2 <- ggplot() +
geom_spatraster(data = rast_deltas_future2) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Future predictions of Total Absolute Cover, \n compared to contemporary predictions (future - contemprary) \n- climate model: IPSL-CM5A-MR (France)")) +
scale_fill_gradient2(low = "orange",
mid = "white" ,
high = "purple" ,
midpoint = 0, limits = c(-2,2), na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist_future2 <- ggplot() +
geom_density(aes(predsFuture2[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of relative cover of ",x,
" \n in future scenario (model = IPSL-CM5A-MR (France))")) +
theme_classic()
} else {
# make a map of predictions w/ current data
resp_map <- ggplot() +
geom_spatraster(data = resp_rast_2) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Predictions of relative cover of ",x,",
\n when all functional groups sum to 1")) +
scale_fill_gradient2(low = "brown",
mid = "wheat" ,
high = "darkgreen" ,
midpoint = 0, limits = c(0,1),
na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist <- ggplot() +
geom_density(aes(contempPreds[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of relative cover of ",x)) +
theme_classic()
# make a map of future preds w/ model 1
resp_map_future1 <- ggplot() +
geom_spatraster(data = rast_deltas_future1) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Future predictions of relative cover of ",x,", \n compared to contemporary predictions (future - contemprary) \n - climate model: BNU-ESM")) +
scale_fill_gradient2(low = "orange",
mid = "white" ,
high = "purple" ,
midpoint = 0, limits = c(-1,1), na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist_future1 <- ggplot() +
geom_density(aes(predsFuture1[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of relative cover of ",x,
" \n in future scenario (model = BNU-ESM)")) +
theme_classic()
# make a map of future preds w/ model 2
resp_map_future2 <- ggplot() +
geom_spatraster(data = rast_deltas_future2) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Future predictions of relative cover of ",
x,
", \n compared to contemporary predictions (future - contemprary) \n- climate model: IPSL-CM5A-MR (France)")) +
scale_fill_gradient2(low = "orange",
mid = "white" ,
high = "purple" ,
midpoint = 0, limits = c(-1,1), na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist_future2 <- ggplot() +
geom_density(aes(predsFuture2[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of relative cover of ",x,
" \n in future scenario (model = IPSL-CM5A-MR (France))")) +
theme_classic()
}
return(list("rast" = resp_rast_2,
"map" = resp_map,
"hist" = resp_hist,
"rast_future1" = resp_rast_future1_2,
"map_future1" = resp_map_future1,
"hist_future1" = resp_hist_future1,
"rast_future2" = resp_rast_future2_2,
"map_future2" = resp_map_future2,
"hist_future2" = resp_hist_future2))
})
names(relCover_A_Maps_contemp) <- c("relCoverA_totalTree", "relCoverA_totalHerb", "relCoverA_shrub", "relCoverA_bareGround",
"relCoverA_NLTree", "relCoverA_BLTree", "relCoverA_C3Grass", "relCoverA_C4Grass", "relCoverA_Forb",
"totalAbsCover_CONUS")ggarrange(
ggarrange(relCover_A_Maps_contemp$relCoverA_totalTree$map,
relCover_A_Maps_contemp$relCoverA_totalTree$map_future1,
relCover_A_Maps_contemp$relCoverA_totalTree$map_future2,
relCover_A_Maps_contemp$relCoverA_totalTree$hist,
relCover_A_Maps_contemp$relCoverA_totalTree$hist_future1, relCover_A_Maps_contemp$relCoverA_totalTree$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_A_Maps_contemp$relCoverA_totalHerb$map,
relCover_A_Maps_contemp$relCoverA_totalHerb$map_future1,
relCover_A_Maps_contemp$relCoverA_totalHerb$map_future2,
relCover_A_Maps_contemp$relCoverA_totalHerb$hist,
relCover_A_Maps_contemp$relCoverA_totalHerb$hist_future1,
relCover_A_Maps_contemp$relCoverA_totalHerb$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_A_Maps_contemp$relCoverA_shrub$map,
relCover_A_Maps_contemp$relCoverA_shrub$map_future1,
relCover_A_Maps_contemp$relCoverA_shrub$map_future2,
relCover_A_Maps_contemp$relCoverA_shrub$hist,
relCover_A_Maps_contemp$relCoverA_shrub$hist_future1,
relCover_A_Maps_contemp$relCoverA_shrub$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_A_Maps_contemp$relCoverA_bareGround$map,
relCover_A_Maps_contemp$relCoverA_bareGround$map_future1,
relCover_A_Maps_contemp$relCoverA_bareGround$map_future2,
relCover_A_Maps_contemp$relCoverA_bareGround$hist,
relCover_A_Maps_contemp$relCoverA_bareGround$hist_future1,
relCover_A_Maps_contemp$relCoverA_bareGround$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_A_Maps_contemp$relCoverA_NLTree$map,
relCover_A_Maps_contemp$relCoverA_NLTree$map_future1,
relCover_A_Maps_contemp$relCoverA_NLTree$map_future2,
relCover_A_Maps_contemp$relCoverA_NLTree$hist,
relCover_A_Maps_contemp$relCoverA_NLTree$hist_future1,
relCover_A_Maps_contemp$relCoverA_NLTree$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_A_Maps_contemp$relCoverA_BLTree$map,
relCover_A_Maps_contemp$relCoverA_BLTree$map_future1,
relCover_A_Maps_contemp$relCoverA_BLTree$map_future2,
relCover_A_Maps_contemp$relCoverA_BLTree$hist,
relCover_A_Maps_contemp$relCoverA_BLTree$hist_future1,
relCover_A_Maps_contemp$relCoverA_BLTree$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_A_Maps_contemp$relCoverA_C3Grass$map,
relCover_A_Maps_contemp$relCoverA_C3Grass$map_future1,
relCover_A_Maps_contemp$relCoverA_C3Grass$map_future2,
relCover_A_Maps_contemp$relCoverA_C3Grass$hist,
relCover_A_Maps_contemp$relCoverA_C3Grass$hist_future1,
relCover_A_Maps_contemp$relCoverA_C3Grass$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_A_Maps_contemp$relCoverA_C4Grass$map,
relCover_A_Maps_contemp$relCoverA_C4Grass$map_future1,
relCover_A_Maps_contemp$relCoverA_C4Grass$map_future2,
relCover_A_Maps_contemp$relCoverA_C4Grass$hist,
relCover_A_Maps_contemp$relCoverA_C4Grass$hist_future1,
relCover_A_Maps_contemp$relCoverA_C4Grass$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_A_Maps_contemp$relCoverA_Forb$map,
relCover_A_Maps_contemp$relCoverA_Forb$map_future1,
relCover_A_Maps_contemp$relCoverA_Forb$map_future2,
relCover_A_Maps_contemp$relCoverA_Forb$hist,
relCover_A_Maps_contemp$relCoverA_Forb$hist_future1,
relCover_A_Maps_contemp$relCoverA_Forb$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_A_Maps_contemp$totalAbsCover_CONUS$map,
relCover_A_Maps_contemp$totalAbsCover_CONUS$map_future1,
relCover_A_Maps_contemp$totalAbsCover_CONUS$map_future2,
relCover_A_Maps_contemp$totalAbsCover_CONUS$hist,
relCover_A_Maps_contemp$totalAbsCover_CONUS$hist_future1,
relCover_A_Maps_contemp$totalAbsCover_CONUS$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ncol = 1
)
Now, show how the relativized cover values change over time
# calculate relativized cover (version w/ all functional groups summing to 1)
predsOverTime <- predsOverTime %>%
mutate(totalAbsCover_CONUS = (absTotalHerb_CONUS + absTotalTree_CONUS + absBareGround_CONUS + absShrub_CONUS)) %>%
mutate(relCoverA_totalTree = absTotalTree_CONUS/totalAbsCover_CONUS,
relCoverA_totalHerb = absTotalHerb_CONUS/totalAbsCover_CONUS,
relCoverA_shrub = absShrub_CONUS/totalAbsCover_CONUS,
relCoverA_bareGround = absBareGround_CONUS/totalAbsCover_CONUS,
relCoverA_NLTree = absNLTree_CONUS/totalAbsCover_CONUS,
relCoverA_BLTree = absBLTree_CONUS/totalAbsCover_CONUS,
relCoverA_C3Grass = absC3grass_CONUS/totalAbsCover_CONUS,
relCoverA_C4Grass = absC4grass_CONUS/totalAbsCover_CONUS,
relCoverA_Forb = absForb_CONUS/totalAbsCover_CONUS
)
# make into a long format
predsOverTime_long <- predsOverTime %>%
select(year,x, y, spatialID, uniqueID, totalAbsCover_CONUS:relCoverA_Forb) %>%
pivot_longer(cols = c(totalAbsCover_CONUS:relCoverA_Forb)) %>%
mutate(year = as.integer(year))
# calculate the coefficient of variation for each cover type at each location
predsOverTime_cv <- predsOverTime_long %>%
group_by(uniqueID, name) %>%
summarize(cv = sd(value)/mean(value))
## add x,y
predsOverTime_cv <- predsOverTime_cv %>%
left_join(predsOverTime %>% select(uniqueID, x, y))
# visualize change over time for a selection of points
predsOverTime_long %>%
filter(uniqueID %in% 1:50) %>%
filter(name %in% c( #"totalAbsCover_CONUS"#,
#"relCoverA_totalTree", "relCoverA_totalHerb",
"relCoverA_shrub" , "relCoverA_bareGround" , "relCoverA_NLTree", "relCoverA_BLTree",
"relCoverA_C3Grass" , "relCoverA_C4Grass", "relCoverA_Forb"
)) %>%
ggplot() +
facet_wrap(~uniqueID) +
geom_area(aes(x = year, y = value, fill = name)) +
# geom_line(data = predsOverTime_long %>%
# filter(uniqueID %in% 1:50) %>%
# filter(name %in% c("relCoverA_totalHerb", "relCoverA_totalTree")),
# aes(x = year, y = value, linetype = name), col = "black") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90)) +
labs(title = "Change in modeled relative cover (all types sum to 1) by functional type from \n 2010 to 2020 for a subset of locations")
predsOverTime_cv %>%
ggplot() +
geom_boxplot(aes(x = name, y = cv, col = name)) +
theme_minimal() +
labs(title = "Coefficients of variation of model-predicted relative cover (all types sum to 1) \n at a given location for each functional group from 2010-2020",
x = "functional group",
y = "coefficient of variation") +
theme(axis.text.x = element_text(angle = 90))
# show the locations of each point that is being used for this temporal analysis
#cropped_states <- cropped_states %>%
# rename(geometry = Shape)
ggplot() +
facet_wrap(~name) +
geom_sf(data = cropped_states, aes()) +
stat_summary2d(data = predsOverTime_cv, aes(x = x, y = y, z = cv), fun = ~ mean(.x, na.rm = TRUE)) +
geom_hex(data = predsOverTime_cv, aes(x = x, y = y, col = cv), stat = "mean", na.rm = TRUE) +
guides(fill = guide_legend(title = "CV")) +
theme_minimal() +
labs(title = "Locations of points included in the random \n sample for examining temporal variation \n color-coded by the coefficient of variation")
Second option: relativize cover so that cover for all groups except trees sums to 1
## for contemporary data
contempPreds <- contempPreds %>%
mutate(total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + #absTotalTree_CONUS +
absBareGround_CONUS + absShrub_CONUS)) %>%
mutate(relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_NLTree = absNLTree_CONUS,
relCoverB_BLTree = absBLTree_CONUS,
relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS
)
## for first climate model
predsFuture1 <- predsFuture1 %>%
mutate(total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + #absTotalTree_CONUS +
absBareGround_CONUS + absShrub_CONUS)) %>%
mutate(relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_NLTree = absNLTree_CONUS,
relCoverB_BLTree = absBLTree_CONUS,
relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS
)
## for second climate model
predsFuture2 <- predsFuture2 %>%
mutate(total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + #absTotalTree_CONUS +
absBareGround_CONUS + absShrub_CONUS)) %>%
mutate(relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_NLTree = absNLTree_CONUS,
relCoverB_BLTree = absBLTree_CONUS,
relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS
)Now, show the results in map form
# list of absolute cover names
responseNames <- c("relCoverB_totalTree", "relCoverB_totalHerb", "relCoverB_shrub", "relCoverB_bareGround",
"relCoverB_NLTree", "relCoverB_BLTree", "relCoverB_C3Grass", "relCoverB_C4Grass", "relCoverB_Forb",
"total_NonTree_AbsCover_CONUS")
relCover_B_Maps_contemp <- lapply(responseNames,
FUN = function(x) {
## contemporary absolute cover
# rasterize response
resp_rast <- contempPreds %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = x,
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
resp_rast_2 <- resp_rast %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
## future model 1 absolute cover
# rasterize response
resp_rast_future1 <- predsFuture1 %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = x,
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
resp_rast_future1_2 <- resp_rast_future1 %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# calculate differences between contemporary preds and future 1 preds (deltas)
rast_deltas_future1 <- resp_rast_future1_2 - resp_rast_2
## future model 2 absolute cover
# rasterize response
resp_rast_future2 <- predsFuture2 %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = x,
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
resp_rast_future2_2 <- resp_rast_future2 %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# calculate differences between contemporary preds and future 1 preds (deltas)
rast_deltas_future2 <- resp_rast_future2_2 - resp_rast_2
if (x == "total_NonTree_AbsCover_CONUS") {
# make a map of predictions w/ current data
resp_map <- ggplot() +
geom_spatraster(data = resp_rast_2) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Predictions of Total Absolute Cover (non-tree) with contemporary data")) +
scale_fill_gradient2(low = "lightblue",
#mid = "wheat" ,
high = "darkblue" ,
#midpoint = 0,
limits = c(0,max(terra::values(resp_rast_2$mean), na.rm= TRUE)),
na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist <- ggplot() +
geom_density(aes(contempPreds[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of Total Absolute Cover (non-tree) with contemporary data")) +
theme_classic()
# make a map of future preds w/ model 1
resp_map_future1 <- ggplot() +
geom_spatraster(data = rast_deltas_future1) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Future predictions of Total Absolute Cover (non-tree) \n compared to contemporary predictions (future - contemprary) \n - climate model: BNU-ESM")) +
scale_fill_gradient2(low = "orange",
mid = "white" ,
high = "purple" ,
midpoint = 0, limits = c(-2,2), na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist_future1 <- ggplot() +
geom_density(aes(predsFuture1[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of relative cover of ",x,
" \n in future scenario (model = BNU-ESM)")) +
theme_classic()
# make a map of future preds w/ model 2
resp_map_future2 <- ggplot() +
geom_spatraster(data = rast_deltas_future2) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Future predictions of Cover (non-tree), \n compared to contemporary predictions (future - contemprary) \n- climate model: IPSL-CM5A-MR (France)")) +
scale_fill_gradient2(low = "orange",
mid = "white" ,
high = "purple" ,
midpoint = 0, limits = c(-2,2), na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist_future2 <- ggplot() +
geom_density(aes(predsFuture2[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of relative cover of ",x,
" \n in future scenario (model = IPSL-CM5A-MR (France))")) +
theme_classic()
} else {
# make a map of predictions w/ current data
resp_map <- ggplot() +
geom_spatraster(data = resp_rast_2) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Predictions of relative cover of ",x,
", \n when all functional groups except trees sum to 1")) +
scale_fill_gradient2(low = "brown",
mid = "wheat" ,
high = "darkgreen" ,
midpoint = 0, limits = c(0,1),
na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist <- ggplot() +
geom_density(aes(contempPreds[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of relative cover of ",x)) +
theme_classic()
# make a map of future preds w/ model 1
resp_map_future1 <- ggplot() +
geom_spatraster(data = rast_deltas_future1) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Future predictions of relative cover of ",x,", \n compared to contemporary predictions (future - contemprary) \n - climate model: BNU-ESM")) +
scale_fill_gradient2(low = "orange",
mid = "white" ,
high = "purple" ,
midpoint = 0, limits = c(-1,1), na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist_future1 <- ggplot() +
geom_density(aes(predsFuture1[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of relative cover of ",x,
" \n in future scenario (model = BNU-ESM)")) +
theme_classic()
# make a map of future preds w/ model 2
resp_map_future2 <- ggplot() +
geom_spatraster(data = rast_deltas_future2) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Future predictions of relative cover of ",
x,
", \n compared to contemporary predictions (future - contemprary) \n- climate model: IPSL-CM5A-MR (France)")) +
scale_fill_gradient2(low = "orange",
mid = "white" ,
high = "purple" ,
midpoint = 0, limits = c(-1,1), na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist_future2 <- ggplot() +
geom_density(aes(predsFuture2[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of relative cover of ",x,
" \n in future scenario (model = IPSL-CM5A-MR (France))")) +
theme_classic()
}
return(list("rast" = resp_rast_2,
"map" = resp_map,
"hist" = resp_hist,
"rast_future1" = resp_rast_future1_2,
"map_future1" = resp_map_future1,
"hist_future1" = resp_hist_future1,
"rast_future2" = resp_rast_future2_2,
"map_future2" = resp_map_future2,
"hist_future2" = resp_hist_future2))
})
names(relCover_B_Maps_contemp) <- c("relCoverB_totalTree", "relCoverB_totalHerb", "relCoverB_shrub", "relCoverB_bareGround",
"relCoverB_NLTree", "relCoverB_BLTree", "relCoverB_C3Grass", "relCoverB_C4Grass", "relCoverB_Forb",
"total_NonTree_AbsCover_CONUS")ggarrange(
ggarrange(relCover_B_Maps_contemp$relCoverB_totalTree$map,
relCover_B_Maps_contemp$relCoverB_totalTree$map_future1,
relCover_B_Maps_contemp$relCoverB_totalTree$map_future2,
relCover_B_Maps_contemp$relCoverB_totalTree$hist,
relCover_B_Maps_contemp$relCoverB_totalTree$hist_future1,
relCover_B_Maps_contemp$relCoverB_totalTree$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_B_Maps_contemp$relCoverB_totalHerb$map,
relCover_B_Maps_contemp$relCoverB_totalHerb$map_future1,
relCover_B_Maps_contemp$relCoverB_totalHerb$map_future2,
relCover_B_Maps_contemp$relCoverB_totalHerb$hist,
relCover_B_Maps_contemp$relCoverB_totalHerb$hist_future1,
relCover_B_Maps_contemp$relCoverB_totalHerb$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_B_Maps_contemp$relCoverB_shrub$map,
relCover_B_Maps_contemp$relCoverB_shrub$map_future1,
relCover_B_Maps_contemp$relCoverB_shrub$map_future2,
relCover_B_Maps_contemp$relCoverB_shrub$hist,
relCover_B_Maps_contemp$relCoverB_shrub$hist_future1,
relCover_B_Maps_contemp$relCoverB_shrub$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_B_Maps_contemp$relCoverB_bareGround$map,
relCover_B_Maps_contemp$relCoverB_bareGround$map_future1,
relCover_B_Maps_contemp$relCoverB_bareGround$map_future2,
relCover_B_Maps_contemp$relCoverB_bareGround$hist,
relCover_B_Maps_contemp$relCoverB_bareGround$hist_future1,
relCover_B_Maps_contemp$relCoverB_bareGround$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_B_Maps_contemp$relCoverB_NLTree$map,
relCover_B_Maps_contemp$relCoverB_NLTree$map_future1,
relCover_B_Maps_contemp$relCoverB_NLTree$map_future2,
relCover_B_Maps_contemp$relCoverB_NLTree$hist,
relCover_B_Maps_contemp$relCoverB_NLTree$hist_future1,
relCover_B_Maps_contemp$relCoverB_NLTree$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_B_Maps_contemp$relCoverB_BLTree$map,
relCover_B_Maps_contemp$relCoverB_BLTree$map_future1,
relCover_B_Maps_contemp$relCoverB_BLTree$map_future2,
relCover_B_Maps_contemp$relCoverB_BLTree$hist,
relCover_B_Maps_contemp$relCoverB_BLTree$hist_future1,
relCover_B_Maps_contemp$relCoverB_BLTree$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_B_Maps_contemp$relCoverB_C3Grass$map,
relCover_B_Maps_contemp$relCoverB_C3Grass$map_future1,
relCover_B_Maps_contemp$relCoverB_C3Grass$map_future2,
relCover_B_Maps_contemp$relCoverB_C3Grass$hist,
relCover_B_Maps_contemp$relCoverB_C3Grass$hist_future1,
relCover_B_Maps_contemp$relCoverB_C3Grass$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_B_Maps_contemp$relCoverB_C4Grass$map,
relCover_B_Maps_contemp$relCoverB_C4Grass$map_future1,
relCover_B_Maps_contemp$relCoverB_C4Grass$map_future2,
relCover_B_Maps_contemp$relCoverB_C4Grass$hist,
relCover_B_Maps_contemp$relCoverB_C4Grass$hist_future1,
relCover_B_Maps_contemp$relCoverB_C4Grass$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_B_Maps_contemp$relCoverB_Forb$map,
relCover_B_Maps_contemp$relCoverB_Forb$map_future1,
relCover_B_Maps_contemp$relCoverB_Forb$map_future2,
relCover_B_Maps_contemp$relCoverB_Forb$hist,
relCover_B_Maps_contemp$relCoverB_Forb$hist_future1,
relCover_B_Maps_contemp$relCoverB_Forb$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$map,
relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$map_future1,
relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$map_future2,
relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$hist,
relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$hist_future1,
relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$hist_future2,
ncol = 3,
nrow = 2,
heights = c(1,.35)
),
ncol = 1
)
Now, show how the relativized cover values change over time
# calculate relativized cover (version w/ all functional groups summing to 1)
predsOverTime <- predsOverTime %>%
mutate(total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + #absTotalTree_CONUS +
absBareGround_CONUS + absShrub_CONUS)) %>%
mutate(relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_NLTree = absNLTree_CONUS,
relCoverB_BLTree = absBLTree_CONUS,
relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS
)
# make into a long format
predsOverTime_long <- predsOverTime %>%
select(year,x, y, spatialID, uniqueID, total_NonTree_AbsCover_CONUS:relCoverB_Forb) %>%
pivot_longer(cols = c(total_NonTree_AbsCover_CONUS:relCoverB_Forb)) %>%
mutate(year = as.integer(year))
# calculate the coefficient of variation for each cover type at each location
predsOverTime_cv <- predsOverTime_long %>%
group_by(uniqueID, name) %>%
summarize(cv = sd(value)/mean(value))
## add x,y
predsOverTime_cv <- predsOverTime_cv %>%
left_join(predsOverTime %>% select(uniqueID, x, y))
# visualize change over time for a selection of points
predsOverTime_long %>%
filter(uniqueID %in% 1:50) %>%
filter(name %in% c( #"total_NonTree_AbsCover_CONUS", "relCoverB_totalTree",
#"relCoverB_totalHerb",
"relCoverB_shrub" , "relCoverB_bareGround" , "relCoverB_NLTree", "relCoverB_BLTree",
"relCoverB_C3Grass" , "relCoverB_C4Grass", "relCoverB_Forb" )) %>%
ggplot() +
facet_wrap(~uniqueID) +
geom_area(aes(x = year, y = value, fill = name)) +
geom_line(data = predsOverTime_long %>%
filter(uniqueID %in% 1:50) %>%
filter(name %in% c("relCoverB_totalHerb", "relCoverB_totalTree")),
aes(x = year, y = value, linetype = name), col = "black") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90)) +
labs(title = "Change in modeled relative cover (all Non-tree types sum to 1) by functional type from \n 2010 to 2020 for a subset of locations")
predsOverTime_cv %>%
ggplot() +
geom_boxplot(aes(x = name, y = cv, col = name)) +
theme_minimal() +
labs(title = "Coefficients of variation of model-predicted relative cover (all Non-tree types sum to 1) \n at a given location for each functional group from 2010-2020",
x = "functional group",
y = "coefficient of variation") +
theme(axis.text.x = element_text(angle = 90))
# show the locations of each point that is being used for this temporal analysis
#cropped_states <- cropped_states %>%
# rename(geometry = Shape)
ggplot() +
facet_wrap(~name) +
geom_sf(data = cropped_states, aes()) +
stat_summary2d(data = predsOverTime_cv, aes(x = x, y = y, z = cv), fun = ~ mean(.x, na.rm = TRUE)) +
geom_hex(data = predsOverTime_cv, aes(x = x, y = y, col = cv), stat = "mean", na.rm = TRUE) +
guides(fill = guide_legend(title = "CV")) +
theme_minimal() +
labs(title = "Locations of points included in the random \n sample for examining temporal variation \n color-coded by the coefficient of variation")